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14.  ABSTRACT 

Further  computational  studies  and  improved  methods  leading  to  prediction  of  properties  and  behavior  of  high-nitrogen  materials 
continue  as  both  computational  capabilities  and  methods  advance.  As  detailed  in  the  2007  book  chapter  titled  “Computational 
Aspects  of  Nitrogen-Rich  HEDMs”  [Rice  BM,  Byrd  EFC,  Mattson  WD.  Struct  Bond.  2007;125:153-194],  a  multitude  of 
models  were  developed  to  predict  key  performance  properties  of  energetic  materials  that  include  numerous  high-nitrogen 
systems  based  on  atomistic  quantum  mechanical  simulations.  In  this  report  we  will  provide  updates  to  these  developments. 
Additionally,  a  number  of  theoretical  studies  have  since  been  undertaken  to  explore  the  structure  and  behavior  of  novel 
condensed  phases  of  polymeric  nitrogen.  This  work  shall  outline  advances  in  the  field  since  its  first  publication  in  2007. 
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1.  Introduction 


A  high-energy-density  material  (HEDM)  must  meet  multiple  various  performance 
and  material  properties  before  selection  for  inclusion  in  a  munition.  These 
properties  range  from  critical  perfonnance  properties,  such  as  heats  of  formation 
and  crystalline  densities  (key  for  determining  detonation  and  propellant 
perfonnance),  to  melting  points,  compatibility  with  other  fonnulation  materials, 
sensitivities  to  insult,  and  environmental  impact.  While  most  research  has  focused 
on  heats  of  formation  and  densities,  researchers  have  begun  investigating  whether 
models  can  be  developed  for  the  prediction  of  other  properties.  In  this  report,  we 
will  present  successes  and  failures  for  several  properties. 

2.  Crystal  Densities 

Beginning  with  crystalline  densities,  there  have  been  numerous  approaches  to 
predict  the  density  without  any  experimental  knowledge  of  the  system  under  study. 
These  methods  can  be  broken  into  2  main  camps:  those  based  on  building  full 
3 -dimensional  (3-D)  representations  of  the  crystal  (ab  initio  crystal  prediction)  and 
those  based  on  models  using  only  the  infonnation  obtained  from  a  single  isolated 
molecule  (molecular  volumes).  Both  of  these  methods  have  been  summarized  in 
Rice  et  al.1  and  both  have  demonstrated  some  success. 

While  ab  initio  crystal  prediction  methods  can  provide  both  predicted  densities  and 
X-ray  spectra,  these  methods  are  usually  computationally  intensive  and/or  typically 
require  the  use  of  classical  force  fields,  which  must  be  constructed  and/or  tested  for 
the  molecular  system  of  interest.  For  a  more  recent  series  of  review  articles  on 
crystal  structure  prediction  methods,  see  Atahan-Evrenk  and  Aspuru-Guzik.2  As 
only  the  density  is  required  for  use  in  predicting  performance  properties,  we  shall 
focus  on  methods  yielding  solely  the  density.  In  our  earlier  work,3  we  initially 
correlated  the  crystal  density  to  a  simple  functional: 


where  M  was  the  molecular  mass  of  the  single  molecule  and  Vm  was  the  volume 
inside  the  0.001  au  isosurface  of  electron  density  surrounding  the  molecule,  which 
was  calculated  using  the  B3LYP/  6-3 1  G**4-6  density  functional  theory  (DFT)  as 
implemented  in  the  Gaussian  program  package.7  Applying  this  method  to  180 
neutral  carbon-hydrogen-nitrogen-oxygen  (CHNO)-containing  molecules  and  38 
high-nitrogen  systems  yielded  root-mean-square  (rms)  percent  deviations  of  3.7% 
and  3.4%,  respectively,  compared  to  experiment.  When  applied  to  ionic  molecules, 
the  resultant  predictions  were  worse,  with  a  6.6%  rms  disagreement  from 
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experiment.  While  applying  a  correction  based  on  the  number  of  hydrogens  in  the 
moiety  improved  the  results  somewhat  (4.8%  rms  error),  later  refinements  to  the 
model,  founded  on  the  works  of  Politzer  et  al.8,9  suggesting  the  addition  of 
corrections  based  on  electrostatic  interactions,  would  provide  an  improved 
correlation.  In  reparametrizing  equations  from  the  Politzer  et  al.  approach,  we 
significantly  improved  the  crystal  density  predictor  with  new  rms  errors  of  2.7% 
for  neutral  molecules  (on  38  test  molecules)  and  3.7%  for  ionic  systems  (on  48  test 
molecules).10  Equations  2  and  3  were  used  for  the  neutral  and  ionic  density 
predictors. 


p  =  «i  (^)  +  A\  +  Yi- 

(2) 

<^(g)+xE(|)+s. 

(3) 

For  the  neutral  molecule  equation,  crt20t  is  the  total  variance  of  the  electrostatic 

potential  mapped  onto  the  0.001  au  isosurface  of  electron  density  of  the  isolated 

molecule;  v  quantifies  the  degree  of  balance  between  the  positive  and  negative 

— + 

potentials  on  the  molecular  surface.  For  the  ionic  molecule  equation,  Fs  is  the 
average  of  the  positive  values  of  the  electrostatic  potential  on  the  0.001  au 
molecular  surface  Vs,  and  A$  is  the  portion  of  the  cation’s  surface  that  has  a  positive 
electrostatic  potential.  Fikewise,  Vs  is  the  average  of  the  negative  values  of  Vs,  and 
A$  is  the  portion  of  the  anion’s  surface  that  has  a  negative  electrostatic  potential. 
These  ratios  are  summed  for  each  charge  moiety  in  the  total  ionic  system.  The  terms 
ai,  Pi,  yi,  a,  P,  y,  and  5  are  all  fitted  constants. 

3.  Solid  Phase  Heats  of  Formation 


The  second  of  the  2  key  performance  properties,  the  solid  phase  heat  of  formation, 
traditionally  has  been  a  computational  challenge.  While  there  exist  numerous 
methods  of  varying  accuracy  and  computation  cost  for  predicting  the  gas  phase 
heats  of  formation,  the  lack  of  an  adequate  treatment  of  the  cohesive  forces  binding 
the  disparate  molecules  together  into  a  solid  has  hampered  the  prediction  of  this 
property.  The  methods  existent  to  date  have  relied  mostly  on  fitted  quantitative 
structure  property  relationships  (QSPRs),  correlating  electronic  structure  properties 
to  experimental  data.  In  Rice  et  al.1,  we  illustrated  a  QSPR  method  for  neutral 
molecules  that  displayed  reasonable  accuracy.  However,  we  demonstrated  there 
were  no  accurate  predictive  tools  for  ionic  compounds.  In  2009,  Byrd  and  Rice11 
examined  a  host  of  different  methods  to  compute  the  gas  phase  heats  of  formation 
along  with  different  models  to  determine  the  lattice  enthalpy,  which,  when  added, 
will  produce  the  desired  solid  phase  heat  of  formation.  Two  methods  were  used  to 
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compute  the  gas  phase  heats  of  formation,  the  very  accurate  (and  computationally 
expensive)  G3MP2B3  theory12,13  and  a  more  empirical  scheme,  while  6  methods 
were  used  to  detennine  the  lattice  enthalpies.  These  methods  ranged  from  extended 
lattice  summations  (requiring  the  crystal  structure,  a  necessary  property  obtained 
through  computationally  costly  crystal  structure  prediction  calculations  for  notional 
compounds)  to  QSPR  methods. 

The  most  recognized  QSPR  model  to  calculate  lattice  enthalpies  used  was  designed 
by  Jenkins  et  al. 14-16  and  was  fitted  to  large  molecular  anions  but  composed  mostly 
of  small  alkali  metal  and  alkaline  earth  cations  (Li,  Na,  K,  Rb,  Cs,  Mg,  Ca,  Sr,  Ba). 
This  restricted  parameterization  is  of  concern  when  applied  to  ionic  high-nitrogen 
materials,  as  most  of  the  anions  are  not  monoatomic  species.  Fortunately,  there 
exists  a  reparametrized  version  of  the  Jenkins  method  for  the  1 : 1  salts,  developed 
by  Gutowski  et  al.,17  that  uses  the  same  functional  fonn  but  corrects  this  deficiency 
in  the  original  fitting  set.  Both  of  these  models  have  an  inverse  cube  root 
dependence  on  the  formula  unit  volume,  which  can  be  computed  as  described 
above. 

The  method  detennined  to  yield  the  most  accurate  results  for  the  condensed  phase 
heat  of  formation  was  the  combination  of  G3MP2B3  theory  (gas  phase  heat  of 
formation)  with  the  Gutowski  method  of  calculating  the  lattice  energy  (24  kcal/mol 
nns  error),  but  this  method  is  restricted  to  1:1  salts.  For  all  other  salts,  the 
G3MP2B3  theory  with  the  Jenkins  method  provided  results  that  were  within 
5  kcal/mol  nns  error  of  the  more  accurate  lattice  summation  methods 
(36.6  kcal/mol  vs.  31.2  kcal/mol  rms  errors).  This  is  fortunate,  as  the  Jenkins  and 
Gutowski  models  require  no  experimental  information,  can  be  rapidly  obtained 
with  very  modest  computational  resources,  and  produce  results  that  are 
approximately  as  accurate  as  those  obtained  using  more  computationally  costly 
methods  that  explicitly  calculate  interatomic  interactions  in  an  ionic  crystal. 

Byrd  and  Rice1 1  also  showed  the  effect  on  predicted  lattice  enthalpies  using  formula 
unit  volumes  based  on  the  refined  densities  using  electrostatic  potential 
information.  Using  the  updated  crystal  density  predictors,  and  then  converting  these 
to  fonnula  unit  volumes  (by  dividing  the  formula  unit  mass  by  the  updated  predicted 
density),  we  observed  minimal  change  in  the  predicted  lattice  enthalpies,  on  the  order 
of  1.4  kcal/mol  average  absolute  difference.9  This  was  expected;  overall,  the  volume 
changes  were  minor,  and  once  filtered  through  an  inverse  cube  root  (the  functional 
fonn  of  the  Jenkins  method),  the  resultant  differences  were  minute. 
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4.  Prediction  of  Vulnerability,  Hazard,  and  Other  Properties 


While  there  has  been  significant  focus  on  the  prediction  of  perfonnance  properties 
of  high-nitrogen  materials,  these  are  not  the  only  factors  to  consider  when  pursuing 
novel  energetic  materials.  Other  concerns  of  note  are  sensitivity  to  insult  (e.g., 
impact,  friction,  electrostatic  discharge),  physical  properties  (e.g.,  melting  points), 
environmental  impact,  and  toxicity  as  well  as  compatibility  with  fonnulation 
ingredients.  These  additional  properties,  however,  usually  involve  complex, 
dynamic  processes  that  do  not  lend  themselves  well  to  “simple”  predictive 
techniques  as  illustrated  previously.  A  prime  example  of  this  is  the  prediction  of 
impact  sensitivity  of  energetic  materials.  There  exists  scores  of  proposed  methods 
to  predict  impact  sensitivities  through  the  correlation  of  molecular  or  material 
properties,  yet  to  date,  none  has  demonstrated  stellar  success.  Our  own  efforts  in 
this  field  have  ranged  from  semi-empirical  QSPR-based  methods  (8  descriptor 
multivariate  linear  regression  equations  yielding  an  R2  of  0.75)18  to  attempted 
correlations  of  electrostatic  properties  on  the  partitioned  quantum  mechanical 
electron  density  derived  from  Bader’s  atoms-in-molecules  technique.19  20  For  these 
studies,  while  many  correlations  were  attempted,  and  several  approximate 
correlations  were  detennined,  all  exhibited  limited  predictive  capability,  and  no 
accurate  quantitative  correlations  were  found.  To  date,  there  exists  no  accurate, 
universal  predictive  methodology  for  sensitivities  for  energetic  materials. 

Another  property  that  has  garnished  considerable  attention  is  the  prediction  of  the 
melting  point  of  a  material.  The  most  common  approach  used  is  that  for  QSPRs, 
which  have  demonstrated  some  successes  for  this  challenging  problem,  some  of 
which  were  developed  for  conventional  energetic  materials.21  In  a  recent  paper,22 
Morrill  and  Byrd  use  the  AMI  semi-empirical  quantum  mechanical  method  to 
design  multivariate  linear  relationships  based  on  the  experimental  melting  points 
for  over  100  high-nitrogen  compounds.  The  models  resulting  from  this  study  had 
R2  values  of  0.94  and  0.90  on  the  test  sets  for  the  2  QSPRs  developed  (nns  errors 
of  ±7.0  and  ±8.8  °C,  respectively),  which  rank  among  the  best  melting  point  models 
published  for  energetic  materials. 

5.  Novel  Polynitrogen  Species  and  Novel  High-Pressure  Phases 
of  Nitrogen 

The  relatively  modest  computational  resources  required  to  obtain  basic  molecular 
information  using  DFT  have  resulted  in  a  proliferation  of  computational  studies  of 
isolated  high-nitrogen  molecules  and  will  not  be  reviewed  here.  However,  the 
advances  in  high-performance  computing,  scalable  DFT  codes,  and  crystal 
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structure  prediction  algorithms  have  allowed  theoretical  explorations  of  structure 
and  behavior  of  novel  high-pressure  phases  of  polymeric  nitrogen  since 
“Computational  Aspects  of  Nitrogen-Rich  HEDMs”1  was  published.  One  of  us 
performed  quantum-based  molecular  dynamic  (QMD)  simulations  in  which  a 
shock  wave  was  initiated  into  a  long  3-D  filament  of  cubic  gauche  nitrogen 
(cg-N)23  through  flyer  plate  impact.  The  shock  response  of  the  filament  was  quite 
different  from  that  which  we  observed  in  QMD  simulations  of  the  shocked 
conventional  energetic  material  pentaerythritol  tetranitrate  (PETN).  In  PETN,  we 
observed  heat  release  reactions  directly  behind  the  shock  front,  while  in  the  shocked 
cg-N  simulations,  a  series  of  energy-absorbing  phase  transfonnations  and 
spontaneous  defect  formations  occurred  immediately  behind  the  shock  front.  The 
decomposition  of  the  polymeric  matrix  occurred  at  the  far  edge  of  the  filament 
opposite  to  the  traveling  shock  front;  the  heat  generated  in  this  decomposition  did 
not  contribute  to  driving  the  front. 

A  follow-up  effort  by  Beaudet  and  Mattson24  focused  on  attempts  to  generate 
ambient  pressure  amorphous  polymeric  nitrogen  condensed  phases,  known  to  exist 
at  low  temperatures  up  to  100  K.25  Generation  of  amorphous  structures  using 
standard  simulated  annealing  methodologies  and  DFT  was  challenging,  leading  us 
to  initiate  the  structure  search  using  semi-amorphous  parts  of  the  highly  defected 
shock  region  in  the  cg-N  simulation.  Thermal  annealing  and  optimization  led  to  the 
discovery  of  a  novel  low-pressure  porous  crystalline  form  of  nitrogen.  The  study 
was  also  particularly  useful  in  identifying  defects  as  sources  of  instability  in  the 
crystal,  infonnation  that  could  be  used  to  guide  stabilization  efforts.  Other  unusual 
structures  of  pure  all-nitrogen  crystals  have  been  theoretically  explored,  including 
a  molecular  crystal  of  Ng  molecules,  whose  main  crystal  forces  are  weak  van  der 
Waals  and  electrostatics;  this  crystal  is  predicted  to  be  metastable  at  ambient 
pressure.26  Crystal  structure  predictions  of  all  nitrogen  systems  at  terapascal  (TPa) 
pressures  produced  several  stable  phases,  with  one  surprising  metallic  phase 
consisting  of  N5  and  N2  moieties  and  exhibiting  ionic  features  and  charge  density 
distortions.27  Another  group  perfonned  a  crystal  structure  prediction  search  that 
produced  an  N10  cage-like  diamonoid  structure  stable  at  pressures  above 
260  gigapascals  (GPa).28 

Heteroatomic  mixtures  of  polymeric  nitrogen  have  also  been  explored  to  increase 
stability  and  otherwise  optimize  material  properties.  Polymeric  nitrogen  forms  are 
found  for  compressed  systems  of  various  hydrogen-nitrogen  mixtures,29-31  sodium 
azide,32,33  CO  with  N2, 34,35  and  N2O36.  Polymeric  nitrogen  has  also  been 
encapsulated  within  nanotubes,  both  carbon37  and  silicon  carbide.38  Finally 
nanotube  crystals  of  nitrogen-nitrogen-oxygen  (NNO)  and  nitrogen-phosphorous- 
oxygen  (NPO)  were  found  through  crystal  structure  prediction  methods.39  While 
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the  NNO  extended  solid  is  not  stable  at  ambient  pressure,  the  NPO  crystal  was 
found  to  be  stable  at  zero  pressure  and  has  an  energy  content  86%  higher  than 
l,3,5,7-tetranitro-l,3,5,7-tetra-azacyclo-octane  (HMX). 

6.  Conclusions 


Increased  interest  in  developing  these  novel  high-nitrogen  materials,  particularly 
the  novel  high-pressure  condensed  phases  of  polymeric  forms  of  nitrogen,  will 
continue  to  benefit  from  theoretical  predictions  of  structure  and  properties  of  the 
materials,  particularly  in  providing  guidance  in  design  and  stabilization. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


6-3 1G** 

Pople’s  Gaussian  double-zeta  polarized  basis  set  with  d  polarization 
functions  on  each  of  the  atoms  Li  through  Ca  and  p  polarization 
functions  on  H  and  He 

V 

balance  parameter 

I 

average  electrostatic  surface  potential 

a2 

variance  of  electrostatic  surface  potential 

P 

crystalline  density 

AMI 

semi-empirical  method 

APG 

Aberdeen  Proving  Ground 

ARL 

US  Army  Research  Laboratory 

au 

atomic  unit 

Ba 

barium 

B3LYP 

Becke  3  parameter  exchange  with  Lee-Yang-Parr  correlation  DFT 
functional 

Ca 

calcium 

cg-N 

cubic  gauche  nitrogen 

CHNO 

carbon-hydrogen-nitrogen-oxygen 

Cs 

cesium 

DFT 

density  functional  theory 

G3MP2B3 

G3  variant  that  uses  MP2  instead  of  MP4  for  the  basis  set  extension 
corrections  and  uses  the  B3LYP  structures  and  frequencies  instead 
of  the  Hartree-Fock  values 

GPa 

gigapascal 

HEDM 

high-energy-density  material 

HMX 

1,3,5 ,7-tetranitro- 1 ,3 ,5 ,7-tetra-azacyclo-octane 

K 

potassium 

kcal/mol 

kilocalories  per  mole  (unit  of  energy) 
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Li 

lithium 

M 

molecular  mass 

Mg 

magnesium 

Na 

sodium 

NNO 

nitrogen-nitrogen-oxygen 

NPO 

nitrogen-phosphorous-oxygen 

PETN 

pentaerythritol  tetranitrate 

QMD 

quantum-based  molecular  dynamics 

QSPR 

quantitative  structure  property  relationships 

R2 

coefficient  of  determination 

Rb 

rubidium 

Sr 

strontium 

TPa 

terapascal 

nns 

root  mean  square 

V 

electrostatic  potential 

Vm 

volume  inside  the  0.001  au  isosurface  of  electron  density 
surrounding  the  molecule 
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